Nonlinear electrochemical relaxation around conductors 



Kevin T. Chu 1 - 2 and Martin Z. Bazant 1 
1 Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02139 
Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, NJ 08544 

(Dated: February 5, 2008) 

We analyze the simplest problem of electrochemical relaxation in more than one dimension - the 
response of an uncharged, ideally polarizable metallic sphere (or cylinder) in a symmetric, binary 
electrolyte to a uniform electric field. In order to go beyond the circuit approximation for thin double 
layers, our analysis is based on the Poisson-Nernst-Planck (PNP) equations of dilute solution theory. 
Unlike most previous studies, however, we focus on the nonlinear regime, where the applied voltage 
across the conductor is larger than the thermal voltage. In such strong electric fields, the classical 
model predicts that the double layer adsorbs enough ions to produce bulk concentration gradients 
and surface conduction. Our analysis begins with a general derivation of surface conservation laws 
in the thin double-layer limit, which provide effective boundary conditions on the quasi-neutral 
bulk. We solve the resulting nonlinear partial differential equations numerically for strong fields 
and also perform a time-dependent asymptotic analysis for weaker fields, where bulk diffusion and 
surface conduction arise as first-order corrections. We also derive various dimensionless parameters 
comparing surface to bulk transport processes, which generalize the Bikerman-Dukhin number. Our 
results have basic relevance for double-layer charging dynamics and nonlinear electrokinetics in the 
ubiquitous PNP approximation. 

PACS numbers: 82.45.Gj, 82.45.Jn, 66.10.-x 
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I. INTRODUCTION 

Diffuse-charge dynamics plays an important role in 
the response of electrochemical and biological systems 
subject to time-dependent voltages or electric fields pj. 
The classical example is impedance spectroscopy in elec- 
trochemistry n n s n, but electrochemical relax- 
ation is also being increasingly exploited in colloids 
and microfluidics [6j- For example, alternating elec- 
tric fields have been used to pump or mix liquid elec- 
trolytes [E H H^, HH [H, Hi Q [IE HE HE El, to 
separate or self-assemble colloids near electrodes [IE I2E 
l2lL 122. l23l I24L |25| . a nd to manipulate polarizable par- 
ticleTl^ppQHipi EE HI H3 or biological cells and 
vesicles I33.l33.l33. 



In this paper, we analyze some simple problems exem- 
plifying the nonlinear response of an electrolyte around 
an ideally polarizable object due to diffusion and electro- 
migration. As shown in Figure ^ we consider the ionic 
relaxation around a metallic sphere (a) or cylinder (b) 
subject to a suddenly applied uniform background elec- 
tric field, as in metallic colloids. Equivalently, we con- 
sider a metallic hemi-sphere (c) or half-cylinder (d) on an 
insulating plane, to understand relaxation around metal- 
lic structures on channel walls in micro-electrochemical 
devices. Although we do not consider fluid flow, our 
analysis of nonlinear electrochemical relaxation is an nec- 
essary first step toward understanding associated prob- 
lems of induced-char ge e lectro-osmosis in the same ge- 
ometries [H EE TIE HE EH, and thus it also has 
relevance for the case of AC electro-osmosis at planar 
electrode arrays 0, IE IE EE El EE EE E3| 

In electrochemistry [3E I2E HE HH and colloid sci- 
ence [2E H3 ! it is common to assume that the charged 




FIG. 1: (a) Schematic diagram of metallic colloidal 
sphere in a binary electrolyte subjected to an ap- 
plied electric field, which has the same relaxation 
as a metallic hemisphere on a flat insulting surface, 
shown in (c). We also consider the analogous two- 
dimensional problems of a metallic cylinder (b), or 
a half cylinder on an insulating plane (d). 



double-layer at a metal surface is very thin and thus 
remains in quasi-equilibrium, even during charging dy- 
namics. As a result, for over a century Q], the stan- 
dard model of electrochemical relaxation has been an 
equivalent circuit, where the neutral bulk is represented 
by an Ohmic resistor and the double layer by a surface 
impedance [E IE El 01 , which reduces to a linear capac- 
itor at an ideally polarizable surface. For our model 
problems, this "RC circuit" model was first applied to 
electrochemical relaxation around a sphere by Simonov 
and Shilov |4l| and around a cylinder by Bazant and 
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Squires [3113. Similar RC-circuit analysis has been ap- 
plied extensively to planar electrode arrays in microflu- 
idic devices, following with Ramos et al. Q and Aj- 
dari Q. 

While convenient for mathematical analysis and often 
sufficiently accurate, circuit models neglect the possi- 
bility of bulk concentration gradients, which can arise 
at large applied volta ges [ l| and/or when the surface is 
highly charged |42ll43lT43.l45| . as well as nonuniform sur- 
face transport of ions through the double layer [ifil l47j . 
Dukhin and Shilov l42l l45| and later Hinch, Sherwood, 
Chew, and Sen |43ll44j made significant progress beyond 
the simple circuit model by including bulk diffusion in 
their studies of double layer polarization around highly 
charged spherical particles (of fixed charged density) in 
weak applied fields. One of the main results of their anal- 
ysis is that for weak applied electric fields, bulk concen- 
tration gradients appear as a small correction (on the or- 
der of the applied electric field) to a uniform background 
concentration. In this work, we will lift the weak-field 
restriction and consider the nonlinear response of the 
system to strong applied fields, using the same mathe- 
matical model as in all prior work - the Poisson-Nernst- 
Planck (PNP) equations of dilute solution theory [3^ ■ As 
we shall see, nonlinear response generally involves non- 
negligible bulk diffusion. 

There are also difficulties with the traditional macro- 
scopic picture of the double layer as an infinitely thin sur- 
face impedance, or possibly some more general nonlinear 
circuit element. At the microscopic level, the double layer 
has a more complicated structure with at least two dif- 
ferent regimes: a diffuse layer where the ions move freely 
in solution and a compact surface layer, where ions may 
be condensed in a Stern mono-layer with its own physical 
features (such as surface ca pac itance, surface diffusivity, 
and surface roughness) [351, [3y, [37J • The surface capaci- 
tance may also include the effect of a dielectric coating, 
where ions do not penetrate @- Mathematically, in one 
dimension the capacitor model can be derived as an ef- 
fective boundary condition for the neutral bulk (Nernst- 
Planck) equations by asymptotic analysis of the PNP 
equations in the thin double-layer limit Extending 
this analysis to higher dimensions, however, requires al- 
lowing for tangential "surface conduction" through the 
double layer on the conductor for large applied electric 
fields. Here, we derive effective boundary conditions for 
the neutral bulk in the PNP model by following a gen- 
eral mathematical method for surface conservation laws 
at microscopically diffuse interfaces [48) . 

The nonlinear response of electrochemical systems to 
strong applied fields seems to be relatively unexplored. 
To our knowledge, the only prior mathematical study 
of nonlinear relaxation with the PNP model is the re- 
cent work of Bazant, Thornton, and Ajdari on the one- 
dimensional problem of parallel-plate blocking electrodes 
applying a sudden pulsed voltage 1]. (The same analy- 
sis has been recently extended to "modified-PNP" models 
accounting for steric effects in concentrated solutions [I^ i 



which would also be an important extension of our work.) 
For applied voltages in the weakly nonlinear regime, 
which is analogous to weak applied fields in our problems, 
they find that the relaxation of the cell to the steady state 
requires bulk diffusion processes that appear as a small 
correction at 0(e), where the small parameter e is the ra- 
tio of the Debye length to a typical scale of the geometry. 
(In contrast, in Refs. |42j, |43J, |4J, |45j, it appears that 
primarily the the strength of the applied electric field is 
assumed to be small, although the mathematical limit 
is not explicitly defined.) For applied voltages in the 
strongly nonlinear regime, they show that bulk concen- 
tration gradients can no longer be considered small, since 
0(1) concentration variations appear. In both regimes, 
the absorption of neutral salt by the double layer (and 
therefore build up of surface ion density) is the key driv- 
ing force for bulk diffusion. This positive salt adsorption 
in response to an applied voltage, first noted in Ref. 0, 
is op pos ite to the classical "Donnan effect" of salt expul- 
sion [I(J, which occurs if the surface chemically injects 
or removes ions during the initial creation of the double 
layer |5(ij| . 

Bazant, et al. also emphasize the importance of both 
the charging and the diffusion time scales in the evolution 
of the electrochemical systems [J . Because circuit mod- 
els inherently neglect diffusion processes, the only char- 
acteristic dynamic time scale that appears is the so-called 
RC charging time, r c = A^L/D, where Ar> is the Debye 
length, L is the system size, and D is the characteristic 
diffusivity of the ions |70| . However, when concentration 
gradients are present, diffusion processes and dynamics 
at the diffusion time scale, tl = L 2 / D may be important. 
Most theoretical analyses of electrochemical systems only 
consider the dynamics at one of these two dominant time 
scales - effectively decoupling the dynamics at the two 
time scales. This decoupling of the dynamics is natural 
when one considers the wide separation in the time scales 
that govern the evolution of the system: ^> t c . An 
interesting discussion of how the two time scales are cou- 
pled using ideas related to time-dependent asymptotic 
matching is presented in Ref. 1] . 

The paper is organized as follows. We begin in sec- 
tion[n]by carefully considering the thin double-layer limit 
of our model problems, which leads to effective boundary 
conditions for the neutral-bulk equations in section II I II 
The most interesting new boundary conditions are sur- 
face conservation laws, whose physical content we discuss 
in detail in section HVl where we also define dimensionless 
parameters governing the importance of various surface 
transport processes. In section we explore the steady 
response to large applied electric fields in our model prob- 
lems; a notable prediction is the formation of recirculat- 
ing bulk diffusion currents coupled to surface transport 
processes. We then turn to relaxation dynamics in the 
three regimes identified by Bazant et al. ^j} , using simi- 
lar methods, albeit with nontrivial modifications for two 
or three dimensions. We begin with the linear response 
to a weak field in section IVII where we obtain exact solu- 
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tions using transform methods for arbitrary double-layer 
thickness and also consider AC electric fields. Next, in 
section rVIII we use boundary-layer methods in space and 
time to analyze "weakly nonlinear" relaxation in some- 
what larger fields, in the asymptotic limit of thin dou- 
ble layers. Finally, in section IVIIII we comment on the 
challenges of "strongly nonlinear" relaxation, where bulk 
diffusion and surface conduction dominate the dynamics. 
We conclude in section llXl bv discussing limitations of the 
PNP model and directions for future research. 



II. MATHEMATICAL MODEL 
A. PNP Initial-Boundary- Value Problem 

As a model problem, we consider the response of an iso- 
lated, ideally polarizable sphere (or cylinder) subjected 
to a uniform, applied electric field, as shown in Figure ^ 
For simplicity, we focus only on a symmetric, binary elec- 
trolyte where both ionic species have the same diffusivity 
and charge number. In order to study nonlinear effects 
and avoid imposing a time scale, we assume that the uni- 
form electric field is suddenly applied at t = 0. 

As in most (if not all) prior work on electrochemical dy- 
namics, we assume the Poisson-Nernst-Planck equations 
of dilute solution theory j3^ , 



dC+ 
~dT 
dC- 



V • | £>VC_ 



z+eD 
kT 

z + eD 
kT 



(1) 

(2) 
(3) 



where D is the diffusivity, z + is the charge number of the 
positively charged species, e is the charge of a proton, k 
is Boltzmann's constant, T is the absolute temperature, 
and e s is the electric permittivity of the solution. As 
usual, we have used the Nernst-Einstein relation to write 
the mobility in terms of the the diffusivity, b = D/kT. 
It is also useful to define the chemical potentials of the 



ions, 



fi± = kT log C± ± z + e& 



(4) 



from which their fluxes are defined as F± = —bC±VfJ,±. 

At the conductor's surface, we assume the same bound- 
ary conditions as in Ref. 0. We adopt a linear surface- 
capacitance condition on the electrostatic potential |5lj| . 



<9$ 

<P + \ Sir = V 
on 



(5) 



where As is a length characterizing the compact-layer 
surface capacitance (e.g. due to a Stern monolayer or 
a thin dielectric coating), where V is the potential of 
the conductor, which is set either externally or by the 
condition of fixed total charge EE Hi (We will focus 



on the case of zero total charge, where symmetry implies 
V = in our simple geometries.) To focus on charging 
dynamics, we also assume an ideally polarizable surface 
with no- flux boundary conditions: 



dn kT dn 
D dC_ _z+eD r 9£ = 
dn kT ~ dn 



(6) 
(7) 



where the direction of the unit normal is taken to point 
inwards towards the center of the sphere (i.e., outwards 
from the region occupied by the electrolyte solution). 
In the far field, we assume that both the concentration 
and potential profiles tend toward their initial conditions, 
given everywhere by 



C±(R,9,<t>,t = 0) = C 
$(R,6,</>,t = 0) = -E R [ 1 



R? 



(8) 

cos 9. (9) 



where a is the radius of the sphere, C is the bulk con- 
centration far away from the conductor, and E Q is the 
applied electric field. For the case of the cylinder, the 
initial conditions for the concentration profile remains 
the same and the initial electric potential takes the form 



<f>(R,9,t = 0) 



E Q R 1 



— \cos0 (10) 



B. Different Contributions to Ion Transport 

The transport equations an d 10 represent conser- 
vation laws for the ionic species where the flux of each 
species is a combination of diffusion and electromigration. 
Because this paper examines ionic fluxes in detail, let us 
take a moment to fix the notation that we shall use to dis- 
cuss different contributions to charge and mass transport. 
All fluxes will be denoted with the by the variable F (or F 
for scalar components of flux). Superscripts will be used 
to distinguish between the diffusion, (d), and electromi- 
gration, (e), contributions to the flux. Subscripts will 
be used to denote the species or quantity with which the 
flux is associated. Finally, normal and tangential compo- 
nents of a flux will denoted through the use of an extra 
subscript: n for the normal component and t for the tan- 
gential component. As examples, F^ = — 2 +^ P C+V ( I > 
represents the flux of the positively charged species due 
to electromigration and F^l — represents the 

flux of the neutral salt concentration, C = (C+ + CL)/2, 
normal to a surface arising from diffusion. Table [I] pro- 
vides a summary of the various bulk fluxes that shall arise 
in our discussion. (We also abuse notation and use the 
same symbol /j, for chemical potential with dimensions 
and scaled to kT.) 
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TABLE I: Summary of notations for for the various ion fluxes 
and chemical potentials, before and after scaling. 





Dimensional Formula a 


Dimensionless Formula 




kT log C± ± z+e$ 


log C± ± <j> 



F± 



-6C*±V/i± 

-DX7C ± 

z,e~D 



-C±V<I> 



-c±Vp,± 
- (Vc± ± c±V0) 
-Vc± 



? (d) 



Z5VC* + ^j^pV* 
-DVC 



fcT 



pV$ 



(Vc + pV0) 

-Vc 
-pV<?i 



F 



(d) 
p 



DVp + 



+ eZ5 
kT~ 

DVp 



c*v$ 



(Vp + cV<jf>) 
-Vp 
— cVci 



"In these formulae, p is half the charge density (not the total 
charge density). Also, we have abused notation and used the same 
variable p for both the dimensional and dimensionless formulae, p 
in the dimensional formulae is equal to C p in the dimensionless 
formulae. 



Dimensionless Formulation 



To facilitate the analysis of the model problem, it is 
convenient to nondimensionalizc the governing equations 
and boundary conditions. Scaling length by the radius of 
the sphere, a, the time to the diffusion time td — a 2 /D, 
and the electric potential by the thermal voltage divided 
by the cation charge number, kT / z + e, the governing 
equations 0-0 become 



dc + 
dc- 
-e 2 V 2 cf> 



= V- (Vc+ + c + V0) (11) 

= V- (Vc_ -c_V0) (12) 
= (c+-c_)/2 (13) 



where c±, cj>, and t are the dimensionless concentra- 
tions, electric potential and time, respectively, the spa- 
tial derivatives are with respect to the dimensionless po- 
sition, and e is the ratio of the Debye screening length, 
Xd = \l ' il'i k Jc '^ to the radius of the sphere. The bound- 



conditions become 



+ 5e 



dn 



dc + ^ ^ d(p 
dn + dn 

dc- d<p 
dn '~ dn 
c±{r,6,<t>,t = Q) 

<P(r,6,<f>,t = 0) 



E n r 1 — I cos 



(14) 

(15) 

(16) 
(17) 

(18) 



In the far field, the dimensionless concentrations 
approach 1 and the electric potential approaches 
—E r cos 9. Note that in nondimensionalizing the sur- 
face capacitance boundary condition (JSJ , we have chosen 
to introduce a new dimensionless parameter 8 = As/ Ad 
which makes it possible to study the effects of the Stern 
layer capacitance independently from the double layer 
thickness [51). 

Because the charge density and neutral salt concen- 
tration are important for understanding the behavior of 
electrochemical transport at high applied fields, it is of- 
ten useful to formulate the governing equations in terms 
of the average concentration, c = (c+ + c_) /2, and half 
the charge density, p — (c + - 
these definitions (|llf> - H13|) can be rewritten as 



i/2 0, ElEa. Using 



dc 

di 
dp 

di 



V ■ (Vc + pV<j>) 

V ■ (Vp + cV</>) 



(19) 

(20) 
(21) 



Throughout our discussion, we shall alternate between 
this formulation and equations Ijlljl - (|13[) depending on 
the context. The initial and boundary conditions for this 
set of equations are easily derived from lfH|l - JT8J|. Here 
we summarize those boundary conditions that change: 



dc dip 
dn ^ dn 
d<t> 



= 



dp i 

dn dn 

cM,<m = o) 
pM,<m = o) 



1 





(22) 

(23) 

(24) 
(25) 



ary conditions at the surface of the sphere and the initial 



For the remainder of this paper, we shall work primar- 
ily with dimensionless equations. On occasion, we will 
mention the dimensional form of various expressions and 
equations to help make their physical interpretation more 
apparent. 



D. Electroneutral Bulk Equations 

In the context of electrokinetics, it is desirable to 
reduce the complexity of the electrochemical transport 
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problem by replacing the PNP equations with a sim- 
pler set of equations that treats the bulk electrolyte and 
the double lay er as separate entities. Circuit models 
E, El have been used extensively to achieve 
this goal by reducing the transport problem to an elec- 
trostatics problem. However, circuit models make the 
rather stringent assumption that bulk concentrations re- 
main uniform at all times. Unfortunately, at high applied 
electric fields, this assumption is no longer valid because 
concentration gradients become important 0. 

In the present analysis, we consider an alternative sim- 
plification of the PNP equations that allows for bulk con- 
centration variations. Since we are interested in colloidal 
systems where particle diameters are on order of microns, 
e is very small which suggests that we consider the thin 
double layer limit (e — > 0). In this limit, the bulk re- 
mains locally electroneutral, so it is acceptable to replace 
Poisson's equation with the local electroneutrality condi- 
tion [38ll53T|: 

Zid = 0- (26) 

i 

For the case of symmetric, binary electrolytes, (|26[1 leads 
to the electroneutral Nernst-Planck equations: 



= V ■ (cV0) (28) 
p = 0. (29) 

Notice that the common practice of modeling the bulk 
electrolyte as a linear resistor obeying Ohm's law, V 2 4> = 
0, arises from these equations if the concentration profile 
is uniform. 

We emphasize that these equations only describe the 
electrochemical system at the macroscopic level (i.e., in 
the bulk region of the solution). The microscopic struc- 
ture within the double layer, where local electroneutral- 
ity breaks down, is completely neglected. Therefore, any 
physical effects of double layer structure can only be 
incorporated into the mathematical model via effective 
boundary conditions. 



well known result.) Because the general form for the ef- 
fective boundary conditions of locally electroneutral elec- 
trochemical systems has not been extensively discussed 
in the literature, we provide a detailed derivation of these 
boundary conditions and discuss some associated dimcn- 
sionless parameters. 

A. Compact-layer Surface Capacitance 

We begin our discussion of effective boundary con- 
ditions by considering the "Stern boundary condition" 
(|14p. which describes the (linear) capacitance of a pos- 
sible compact layer on the surface. Because Eq. I|14|) 
involves the electric potential and electric field at the in- 
ner edge of the diffuse charge layer, it cannot be directly 
used as a boundary condition for the locally electroneu- 
tral equations l|27|l - l|29|l . However, by rearranging l|I4|) 
and using the GCS model, it is is possible to rewrite the 
Stern boundary condition so that it only explicitly in- 
volves the electric potential and average concentration at 
the macroscopic surface of the electrode (i.e., the outer 
edge of the diffuse charge layer) 0, : 

C + 2<Vcsinh(C/2) =v-<j>. (30) 

Here £ is the potential drop across the diffuse part of the 
double layer (i.e., zeta-potential), v is the potential of 
the metal sphere in the thermal voltage scale, and <f> and 
c are the values of the electric potential and average con- 
centration at the outer edge of the diffuse charge layer. 
With dimensions, l|5U|) is given by 

C + «V^**(^)-v-. (3D 

where £ is t ne dimensional zeta-potential. Note that 
the GCS model for the double layer leads to the conclu- 
sion that the only dependence of the normal derivative 
of the electric potential on the double layer structure is 
through the zeta-potential. A more detailed derivation of 
this form of the Stern boundary condition can be found 
in |51| . 



III. EFFECTIVE BOUNDARY CONDITIONS 
OUTSIDE THE DOUBLE LAYER 

When local electroneutrality is used in place of Pois- 
son's equation, the physical boundary conditions im- 
posed at at electrode surfaces must be modified to ac- 
count for the microscopic structure within the double 
layer. Fortunately, in the e — > limit, the double layer re- 
mains in quasi-equilibrium, so the Gouy-Chapman-Stern 
model can be used to derive effective boundary condi- 
tions for the system. We emphasize that the GCS model 
is not an assumption; rather, it emerges as the leading- 
order approximation in an asymptotic analysis of the thin 
double-layer limit. (See Ref. [jj] for the history of this 



B. Diffuse-layer Surface Conservation Laws 

The effective boundary conditions for ionic fluxes are 
more complicated. Because the physical domain for the 
macroscopic equations ()27(l - Q29JI excludes the diffuse 
part of the double layer, the no-flux boundary conditions 
do not apply - it is possible for there to exist ion flux 
between the bulk region and the double layer. More- 
over, there is also the possibility of ion transport within 
the double layer itself (often neglected) which must be 
accounted for. 

The derivation of the effective flux boundary condi- 
tions (|39f) is based on a general theory of surface conser- 
vation laws at microscopically diffuse interfaces, which 
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FIG. 2: Schematic diagram of normal and tangential 
fluxes involved in the surface conservation law 13211 . 
The shaded circle represents the metal particle; the 
dotted circle represents the outer "edge" of the double 
layer. 



we develop in Ref. |4g. The basic physical idea is to in- 
tegrate out the spatial variation within the double layer 
in the direction normal to the electrode-electrolyte inter- 
face. While intuitively obvious, carrying out the integra- 
tion involves careful use of asymptotic analysis to address 
the mathematical subtleties of integration over boundary 
layers. The theory tells us that effective flux boundary 
conditions have the physically apparent form 



ar 

at 



-*r = -V a -J t ,, 



F„. 



(32) 



where (V s -) denotes surface divergence, 3t,i is the tan- 
gential flux within and F n ^ is the normal flux into the 
boundary layer (see Figure [21 . The effective fluxes J t ,i 
and F Ut i are directly related to the flux P, for the trans- 
port process via 



Jt,i = e 



dz 



(33) 
(34) 



where F and F denote the flux within the boundary layer 
and the flux in the bulk just outside of the boundary layer 
and the integration is over the entire boundary layer (i.e., 
z is the inner variable of a boundary layer analysis). For 
electrochemical transport, the flux is given by the Nernst- 
Planck equation 



(Vc,; + ZiCiV0) 



(35) 



Substituting this expression into l|33|l - l|34|l , rearranging 
a bit, and u sing th e definition of a surface excess concen- 
tration 39, 40,1431 



rv = e 



7idz = e (cj - Ci) dz, 



(36) 



we obtain 

J M = - ^V s r, + ZiTi\J s <p + ezi J Ci\/si>dz)('->7) 



F„ 



dc{ d(j) 

' "o %i c i Ti ■ 

On on 



(38) 



where ip = <fi — (f> is the excess electric potential within the 
boundary layer and V s denotes a surface gradient. As be- 
fore, the tilde (~) and hat (') accents denote the quantities 
within boundary layer and the in the bulk immediately 
outside of the boundary layer. Finally, the effective flux 
boundary conditions for electrochemical transport follow 
by using these results in 1)32(1: 



dTi 
dt 



\7 s T l + ZjYi'SJ s 4> + ezi / Ci\7 s ip dz 
Jo 

dci d<j> 



On On 



(39) 



Note that even though our choice of electric potential 
scale eliminates the need to explicitly refer to the ionic 
charge numbers, 2j, we opt to continue using them in 
the present discussion so that it is clear where the charge 
number should appear for alternative choices of the elec- 
tric potential scale; in the following discussion, the Zi 
are essentially the sign of the "dimensional" ionic charge 
numbers. 

There are a few important features of (|39|l worth men- 
tioning. First, the surface transport term (first term on 
the right hand side) does not always contribute to the 
leading order effective flux boundary condition. Whether 
the surface conduction term must be retained at leading 
order depends on the magnitudes of Ti (which in turn de- 
pends on the zeta-potential) and the tangential compo- 
nent of the bulk electric field. Interestingly, when the the 
surface transport term is significant, the flux boundary 
condition depends explicitly on the small parameter e. 
Also, (|39l) allow for two important physical effects that 
only arise for 2- and 3-D systems: (i) non-uniform dou- 
ble layer charging and (ii) surface transport within the 
double layer itself. The presence of these effects lead to 
richer behavior for 2- and 3-D systems compared to the 
ID system studied in 

To put into a more useful form, we use the GCS 
model of the double layer to express the surface flux den- 
sities in terms of the zeta-potential and the bulk concen- 
tration. From the GCS model HHH], we know that the 
excess concentration of each ionic species is given by 



and that 



7« = c i 



dz 



(e-^ - f) 



sinh 



Z+-0 



(40) 



(41) 



Using these expressions, it is straightforward to show that 
the surface excess concentration is 



2ey 



-ZiC/2 



1 



(42) 
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Therefore, the first two surface flux density terms in 1391) 
can be written as 

r, 

V s Ti+ ZiTiVgC/) = yV s logC + ZiTiVsCj) 



n(z t )Vde- z *V 2 V s (. (43) 



e sg 



To evaluate the last term in the surface flux density, we 
observe that 



dz 



2yfc sinh(C/2) 2 



V s logc 



(44) 



which follows directly by comparing the normal and sur- 
face derivatives of the leading order expression for the 
electric potential within the double layer Q, 0, [S^] 



<ip(z) = 4tanh _1 (tanh(C/4)e 



(45) 



Using this result, the integral in (|39|l greatly simplifies 
and yields 



ezi I Ci\7 s ip dz 
lo 



e sgn(z, 



c e 



yV s logc. 



(46) 



Finally, combining l|43|) and the effective flux 

boundary condition (|39|) becomes 



dt 



= v., 



TiV s logc + ZiTiVgf, 
dci dej) 

— + ZiCi — 

On on 



(47) 



where (dropping hats) c and 4> are understood to be in 
the bulk, just outside the double layer. Notice that the 
tangential gradients in the zeta-potential have vanished 
in this equation so that the surface flux density of the 
individual species is independent of V S C- 

In terms of the (dimcnsionless) chemical potentials of 
the ions, 



/A, = log Cj + Zi 



(48) 



the surface conservation law (|47|l reduces to a very simple 
form, 

dT 

= V s • (IW/Xi) - n ■ CiVfMi (49) 

ot 

where it is clear that the tangential gradient in the bulk 
chemical potential just outside the double layer drives the 
transport of the surface excess concentration of each ion, 
as well as the bulk transport. This form of the surface 
conservation law should hold more generally, even when 
the chemical potential does not come from dilute solution 
theory (PNP equations). 

Effective boundary conditions similar to (|47|) describ- 
ing the dynamics of the double layer have been known 



for some time. In the late 1960s, Dukhin, Deryagin, and 
Shilov essentially used (|47|l in their studies of surface 
conductance and the polarization of the diffuse charge 
layer around spherical par ticles with thin double layers 
at weak applied fields |42j, |45j, |47J . Later, Hinch, Sher- 
wood, Chew, and Sen used similar boundary conditions 
in their extension of the work of Dukhin et al. to explic- 
itly calculate the tangential flux terms for a range of large 
surface potentials and asymmetric electrolytes 0i l^j - 
A key feature of both of these studies is the focus on 
small deviations from bulk equilibrium. As a result, the 
effective boundary conditions used are basically applica- 
tions of (|47() for weak perturbations to the background 
concentration and electric potential. Our work differs 
from these previous analyses because we do not require 
that the bulk concentration only weakly deviates from a 
uniform profile, and we more rigorously justify the ap- 
proximation through the use of matched asymptotics. In 
addition, our analysis is not restricted to the use of the 
GCS model for the double layer; equations l|3T|) - (fH^f) 
are valid even for more general models of the boundary 
layer. 

Before moving on, we write the effective boundary con- 
ditions in dimensional form to emphasize their physical 
interpretation: 



dt 



ZitD 
kT 



IW S $ 



D 



dd 
dn 



Zl eD gg 
kT 1 dn 



(50) 



where I\ is the dimensional surface excess concentration 
of species i and is defined as 



{pi ~ a) 



d dZ 



(51) 



where the integration is only over the double layer. With 
the units replaced, it becomes clear that the effective 
boundary conditions is a two-dimensional conservation 
law for the excess surface concentration Fj with a driv- 
ing force for the flux and a source term that depend only 
on the dynamics away from the surface. Thus, the effec- 
tive boundary conditions naturally generalize the simple 
capacitor picture of the double layer to allow for flow of 
ionic species tangentially along the electrode surface. 



C. Surface Charge and Excess Salt Concentration 

Since the governing equations ljT7ll - l|2~9~|l are formu- 
lated in terms of the salt concentration and charge den- 
sity, it is convenient to derive boundary conditions that 
are directly related to these quantities. Toward this end, 
we define eq and ew to be the surface charge density and 
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surface excess salt concentration, respectively: 

/■OC -j^ /*oo 

q = pdz= - ( 7+ - 7 _)dz 



(52) 



(c — c) <iz 



1 

2 J(I 



( 7+ +7-)«fo. (53) 



By integrating the expressions for the diffuse layer salt 
concentration and charge density 0, 113, 13 , 



c 



c cosh if} 
— csmhip 



(54) 
(55) 



and using l|41|) , both q and w can be expressed as simple 
functions of the zeta-potential and the bulk concentration 
just outside of the double layer: 



q = -2V5sinh(C/2) 
w = 4V£sinh 2 (C/4). 



(56) 
(57) 



Thus, we can combine the effective flux boundary condi- 
tions for individual ions l|39fl to obtain 



dq _ 

:— = ev s 
dt 



V s <? + wV s 



dn 



dw _ 
dt 



V s u> + q\7 s 



cV s ip dz 



pV s ip dz 



dc 
dn 



(58) 



(59) 



Notice that, as in the PNP equations written in terms of 
c and p, there is a symmetry between q and w in these 
equations. As in the previous section, we can use the 
GCS model to rewrite l|58|l and (|59ll solely in terms of 
bulk field variables and the zeta-potcntial: 



dq 

dw 
~dt 



eV s 



(q\7 s log c + w\7 s 
( wS/ s log c + gV s 



4 (60) 
dn 

dc 
dn 



(61) 



which is the form of the effective flux boundary condi- 
tions we shall use in our analysis below. 



IV. SURFACE TRANSPORT PROCESSES 

Before proceeding to the analysis, we discuss the rel- 
ative importance of the various surface transport pro- 
cesses, compared to their neutral bulk counterparts. For 
clarity, we summarize our results for tangential surface 
fluxes in Table ITT1 with and without dimensions. As with 
the associated bulk fluxes summarized in Table [IJ we in- 
troduce different notations for contributions by diffusion 
and electromigration (superscripts (d) and (e), respec- 
tively) to the tangential (subscript t) surface fluxes of 



cations, anions, charge, and neutral salt (subscripts +, 
— , q, and w, respectively). This allows us to define di- 
mensionless parameters comparing the various contribu- 
tions. 

J. J. Bikerman pioneered the experimental and the- 
oretical study of double-layer surface transport 0, |5f| 
and first defined the dimensionless ratio of surface cur- 
rent to bulk current, across a geometrical length scale, a, 
for a uniformly charged double layer and a uniform bulk 
solution |4(j . Using our notation, the Bikerman number 
is 



Bi 



Jt.q 

aJ 



(62) 



where Jo = z+eFo is a reference bulk current in terms of 
the typical diffusive flux, Fq — DCo/a. B. V. Deryagin 
and S. S. Dukhin later added contributions from electro- 
osmotic flow, using the GCS model of the double layer 
(as did Bikerman) (4^, and Dukhin and collaborators 
then used this model to study electrophoresis of highly 
charged particles (with large, but nearly uniform sur- 
face charge) ^2|>0>|5(j. As in other situations, surface- 
conduction effects in electrophoresis are controlled by Bi, 
which thus came to be known in the Russian literature 
as the "Dukhin number", Du. (Dukhin himself denoted 
it by Rel). 

Recently, Bazant, Thornton and Ajdari pointed out 
that the (steady) Bikerman-Dukhin number is equal to 
the ratio of the excess surface salt concentration to its 
bulk counterpart at the geometrical scale a, 



Bi 



r. 



2aCo 



(63) 



and they showed its importance in a one-dimensional 
problem of electrochemical relaxation between parallcl- 
plate electrodes lj. This surprising equivalence (in light 
of the definition of Bi) demonstrates that surface con- 
duction becomes important relative to bulk conduction 
simply because a significant number of ions are adsorbed 
in the double layer compared to the nearby bulk solution. 
This means that salt adsorption leading to bulk diffusion 
should generally occur at the same time as surface con- 
duction, if the double layer becomes highly charged dur- 
ing the response to an applied field or voltage, as in our 
model problems below. 

Unlike prior work, however, our multi-dimensional 
nonlinear analysis allows for nonuniform, time-dependent 
charging of the double layer, so the Bikerman-Dukhin 
number can only be defined locally and in principle could 
vary wildly across the surface. Moreover, since we sepa- 
rate the contributions to surface transport from diffusion 
and electromigration, we can define some new dimen- 
sionless numbers. From Eqs. H6U|I - (|61|I . we find that the 
following surface-to-bulk flux ratios are related to the ex- 
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TABLE II: Summary of surface flux formulae for electrochemical transport. 



Dimensional Formula"' ,c Dimensionless Formula 

h,± — &r±v s ^± -r±v s ^± 

eD , 



- [ot± v s log (c/c ) ± v fl $j - (r± v s log c ± r± v a 

J^ d i -Dr ± v s io g (c/c ) ~ -r±v s io gc 

T (e) T H sfl 



=f^f±v s $ Tr±v s 



Jt,« - [-DQV S log (C/C ) + ^-WV S ®\ -e (gV s log c + ™V S 

jgj -DQV S log (C/C) - eg V s logc 

Jt,™ - [WV S log (C/C ) + i±|^QV s *] -e(™V s logc + qV s 

jg> -DWV S log (C/C ) -«uV s logc 

a We have abused notation and used the same variable T; for both 
the dimensional and dimensionless formulae. r\ in the dimensional 
formulae is equal to C aYi in the dimensionless formulae. 

''The dimensional surface charge density, Q, is defined by Q = 

The dimensional excess neutral salt concentration, W, is defined 
by W = C ^d w - 



cess surface-to-bulk ratio of neutral salt concentration 



= ew = 4- 



h^g I 
aJo 



aF 



(64) 
(65) 



Note that when surface diffusion is neglected, as in most 
prior work, then a = Bi, since surface currents arise 
from electromigration alone. The other surface-to-bulk 
flux ratios are given by the surface-to-bulk ratio of the 
charge density, 







A; 



\4S\ 

aJ 



sinh 



-eC 



l^6,«;l 



(66) 
(67) 



For thin double layers (An < a), we see that the surface- 
to-bulk flux ratios, a and (3, are only significant when £ 
significantly exceeds the thermal voltage kT/e. 

To better understand how the charge and neutral-salt 
fluxes are carried, it is instructive to form the ratio of 
these numbers, 



= tanh 



\Jt,w I 



4kT J 

1 45 1 



i j 



t,w\ 



I J, 



(68) 
(69) 



For weakly charged double layers, z+eC <§C 4/cT, where 
a <§; /3, the (small) surface flux of salt is dominated by 
electromigration, while the surface flux of charge (sur- 
face current) is dominated by diffusion (if bulk concen- 
tration gradients exist). For highly charged double layers 
z + eQ > AkT where a ~ /3, the contributions to each flux 
by diffusion and electromigration become comparable, as 
counterions are completely expelled (q ~ w). 



V. NONLINEAR STEADY RESPONSE FOR 
THIN DOUBLE LAYERS 

A. Effective Equations 

Using the mathematical model developed in the pre- 
vious section, we now examine the steady response of 
a metal sphere or cylinder with thin double layers sub- 
jected to a large, uniform applied electric field. At 
steady-state, the unsteady term is eliminated from the 
governing equations (j2*?|l - > so we nave 



= 
= 



V 2 c 



V • (cV0) . 



(70) 
(71) 



Similarly, the flux boundary conditions l|6(JI) - (|61l) bc- 



t,q 



= eV s ■ (gV s logc + u/V s <; 
= eV s ■ (wV s logc-f-gV s <; 



dn 
dc 

dn 



(72) 
(73) 
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Notice that we have retained the surface transport terms 
in these boundary conditions even though they appear at 
0(e). At large applied fields, it is no longer valid to order 
the terms in an asymptotic expansions by e alone. We 
also need to consider factors of the form ee^ or esinh(C) 
since e 1 * may be as large as 0(1/ eE) for large applied 
fields. Since both q and w contain factors which grow 
exponentially with the zeta-potential, we cannot discard 
the surface conduction terms in and (J73J). Finally, 
the Stern boundary condition 1)30(1 remains unchanged 
because it does not involve any time derivatives. 

As mentioned earlier, the steady problem exhibits in- 
teresting features that have not been extensively ex- 
plored. Unfortunately, the nonlinearities present in the 
governing equations and boundary conditions make it dif- 
ficult to proceed analytically, so we use numerical meth- 
ods to gain insight into the behavior of the system. In 
this section, we first briefly describe the numerical model 
we use to study the system. We then use the numerical 
model to study the development of 0(1) bulk concen- 
tration variations and their impact on transport around 
metal colloid spheres. 



B. Numerical Model 

To solve (|70t - <|71l) numerically in a computationally 
efficient manner, we use a pseudospectral method |57l 
IBsl l59l ]. For problems in electrochemical transport, 
pseudospectral methods are particularly powerful be- 
cause they naturally resolve boundary layers by plac- 
ing more grid points near boundaries of the physical do- 
main l58l . We further reduces the computational 
complexity and cost of the numerical model by taking 
advantage of the axisymmetry of the problem to reduce 
the numerical model to two dimensions (as opposed to 
using a fully 3-D description). 

For the computational grid, we use a tensor product 
grid of a uniformly spaced grid for the polar angle direc- 
tion and a shifted semi-infinite rational Chebyshev grid 
for the radial direction |Hi| . To obtain the discretized 
form of the differential operators on this grid, we use 
Kronecker products of the differentiation matrices for 
the individual one-dimensional grids [57| . The numeri- 
cal model is then easily constructed using collocation by 
replacing field variables and continuous operators in the 
mathematical model by grid functions and discrete oper- 
ators. The resulting nonlinear, algebraic system of equa- 
tions for the values of the unknowns at the grid points 
is solved using a standard Newton iteration with contin- 
uation in the strength of the applied electric field. The 
Jacobian for the Newton iteration is computed exactly 
by using a set of simple matrix-based differentiation rules 
derived in the appendix of (53|. By using the exact Ja- 
cobian, the convergence rate of the Newton iteration is 
kept low; typically less than five iterations are required 
for each value of the continuation parameter before the 
residual of the solution to the discrete system of equa- 



tions is reduced to an absolute tolerance of 10 -8 fn\- 
Directly computing the Jacobian in matrix form had the 
additional advantage of making it easy to implement the 
numerical model in MATLAB, a high-level programming 
language with a large library of built-in functions. It is 
also worth mentioning that we avoid the problem of deal- 
ing with infinite values of the electric potential by formu- 
lating the numerical model in terms of <j) + Er cos 9, the 
deviation of the electric potential from that of the uni- 
form applied electric field, rather than <fi itself. 

In our numerical investigations, we used the numerical 
method described above with 90 radial and 75 angular 
grid points. This grid resolution balanced the combined 
goals of high accuracy and good computational perfor- 
mance. Figure |3| shows typical solutions for the concen- 
tration and electric potential (relative to the background 
applied potential) for large applied electric fields. A com- 
parison of the concentration and electric potential at the 
surface of the sphere for varying values of E are shown 
in Figure 




FIG. 3: Numerical solutions for the concentration c 
(left) and excess electric potential ip (right) for E — 10, 
e = 0.01, 5 = 1. Notice the large gradients in the 
concentration profile near the surface of the sphere. 



C. Enhanced Surface Excess Concentration and 
Surface Conduction 

Perhaps the most important aspects of the steady re- 
sponse at high applied electric fields are the enhanced 




FIG. 4: Bulk concentration c and electric potential (f> 
at the surface of the sphere for varying values of the 
applied electric field. In these figures, e = 0.01 and 

8 = 1. Notice that for E = 15, the poles (9 — and 

9 = 7r) are about to be depleted of ions (i.e., c ~ 0). 



11 




FIG. 5: Surface charge density eq (left) and excess 
surface concentration of neutral salt ew (right) and for 
varying values of the applied electric field. In these 
figures, e = 0.01 and 5 = 1. Notice that for large 
applied fields, ew = 0(1/ E) and eq = 0(1/E) so that 
the surface conduction terms in 1601 and 16H are 0(1). 




e e 

FIG. 6: Tangential surface fluxes for the surface 
charge density |Jt, 9 | (left) and excess surface concen- 
tration of neutral salt |Jt, lu | (right) for varying values 
of the applied electric field. In these figures, e = 0.01 
and 5 = 1. Notice that for large applied fields, the 
surface fluxes are O(l) quantities and make a non- 
negligible leading-order contribution in iltilH and JHB- 



surface excess ion concentration and surface transport 
within the double layer. As shown in Figure |SJ at 
high applied fields, the excess surface concentrations is 
0(l/E), so surface transport within the double layer 
(shown in Figure [BJ becomes non- negligible in the lead- 
ing order effective flux boundary conditions 172|) - (|73|l . 
Interestingly, surface conduction, and are the 

dominant contributions to surface transport (see Fig- 
ure UJ. While there are clearly surface gradients in con- 
centration, surface diffusion is smaller than surface mi- 
gration by a factor on the order of 1/cE. Also, we reit- 
erate that the driving force for surface transport is solely 
from surface gradients of the bulk concentration and bulk 
electric potential; gradients in the zeta-potential do not 
play a role because they are completely canceled out. 

An important feature of the surface fluxes, 3t,q and 
3t,w, shown in Figure El is that they are non-uniform. 
This non-uniformity is strongly influenced by the non- 
uniformity in the tangential electric field (see Figure [SJ. 
The non-uniformity of the surface excess salt concentra- 
tion and charge surface density also play a role but to 
a lesser extent. For the steady problem, the surface ex- 




FIG. 7: Comparison of the mag nitudesof |J (d) | (solid 
lines) and |J' e '| (dashed lines) for the excess surface 
fluxes of eq (left) and ew (right) for an applied electric 
field value of E — 15. In these figures, e = 0.01 and 
5 = 1. Notice that in both cases, the surface diffusion 
is on the order of 1/cE times the surface conduction. 



20 




4 



e 

FIG. 8: Tangential component of bulk electric field, 
\Et\ at surface of sphere for varying values of the ap- 
plied electric field. In these figures, e = 0.01 and 5 = 1. 



cess concentration of ions remains constant in time, so 
the non-uniformity in the surface fluxes leads to non- 
uniform normal fluxes of current and neutral salt from 
the bulk into the double layer (see Figure |5J • Notice 
that the normal flux of neutral salt into the double layer, 
which is given by —dc/dn, shows an injection of neutral 
salt at the poles {—dc/dn > 0), where the charging is 
strongest, and an ejection of neutral salt at the equator 
{—dc/dn < 0), where there is essentially no excess neu- 




FIG. 9: Normal flux of current (left) and neutral 
salt (right) into the double layer for varying values 
of the applied electric field. In these figures, e = 0.01 
and 5 = 1. 
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FIG. 10: Diffusion currents drive transport of neu- 
tral salt near the surface of the sphere. In this figure, 
E — 10, e = 0.01 and 8 — 1. Notice that streamlines 
of neutral salt are closed; current lines start on the 
surface of the sphere near the equator and end closer 
to the poles. 




1 e = it/2 

V 




*. •. e = it/4 













FIG. 11: Comparison of the magnitudes of bulk elec- 
tromigration |Fp e '| (solid lines) and diffusion |Fe | 
(dashed lines) fluxes as a function of distance from 
the surface of the sphere at 6 — 0, n/4, and n/2 (left). 
The figure on the right zooms in on the diffusion flux. 
In these figure, E = 10, e = 0.01 and 5 = 1. No- 
tice that magnitude of the electromigration dominates 
diffusion and that the electromigration term itself be- 
comes dominated by the contribution from the applied 
electric field a short distance from the surface of the 
sphere. 



tral salt build up. This configuration of fluxes leads to 
the neutral salt depletion at the poles and accumulation 
at the equator shown in Figures and ^ Similarly, the 
normal current density —cd(j)/dn, shows an influx of neg- 
ative (positive) current density at the north (south) pole 
and a positive (negative) current density closer to the 
equator. At the equator, there is no normal current den- 
sity because the normal diffusion currents of cations and 
anions exactly balance and there is no normal electric 
field to drive a migration current. 

D. Bulk Diffusion and Concentration Gradients 

One major consequence of surface conduction is trans- 
port of large amounts of neutral salt into the double 
layer. These cause strong concentration gradients to ap- 
pear near the surface of the sphere (Figures |3] and 0}, 
indicating that the usual assumption of a uniform con- 
centration profile is invalid at high electric fields. The 
presence of these large concentration gradients at rela- 
tively low electric fields (E w 5) should not be surpris- 
ing since it is well-known that large voltage effects often 
begin with voltages as low as a few times the thermal 
voltage 1] . The dramatic influence of the voltage arises 
from the exponential dependence of double layer concen- 
trations on the zeta-potential. 

Since the transport of neutral salt is driven by concen- 
tration gradients, the presence of these strong variations 
leads to diffusion currents (see Figure HUIl . An important 
feature of these diffusion currents is that they are closed; 
current lines start on the surface of the sphere near the 
equator where neutral salt is ejected into the bulk (as a 
result of neutral salt transport within the double layer) 
and end close to the poles where neutral salt is absorbed 
by the double layer. These recirculation currents are im- 
portant because they allow the system to conserve the to- 



tal number of cation and anions locally. Without them, 
the local depletion and accumulation of ions would re- 
quire global changes to the bulk concentration (i.e., the 
concentration at infinity would be affected). 

While the presence of diffusion currents is interesting, 
we must be careful in how they are interpreted in terms 
of the motion of individual ion molecules. In actuality, no 
ions are moving purely under the influence of diffusion. 
Rather, the cation and anion migration flux densities are 
slightly imbalanced due to the presence of a concentra- 
tion gradient which results in a net transport of neutral 
salt concentration. 



E. Individual Ion Currents 

Since cations and anions are the physical entities that 
are transported through the electrolyte, it is useful to 
consider the cation and anion flux densities individually. 
As shown in Figure ITT1 in the bulk, the contribution of 
electromigration to transport dominates diffusion. More- 
over, within a short distance from the sphere, the electro- 
migration term itself becomes dominated by the contri- 
bution from the applied electric field. Thus, the concen- 
tration gradient only serves to slightly bias the flux den- 
sities so that cation (anion) motion is slightly retarded 
near the north (south) pole. 

It is more interesting to consider the surface transport 
of the individual ions within the double layer. In the 
northern hemisphere, the double layer is dominated by 
anions; similarly, the southern hemisphere is dominated 
by cations. As a result, transport in each hemisphere is 
primarily due to only one species (see Figure ^J. This 
observation provides a direct explanation for the deple- 
tion and accumulation regions in the concentration pro- 
file in terms of the the motion of individual ions. As 
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FIG. 12: Tangential surface fluxes for the cation | Jt,+ | 
(left) and anion |Jt,-| (right) for varying values of the 
applied electric field. In these figures, e = 0.01 and 
5 = 1. Notice that for large applied fields, the surface 
fluxes are O(l) quantities (in the appropriate hemi- 
sphere) which leads to a non-negligible leading-order 
contribution in 11471 . 



mentioned earlier, the transport is from the poles to the 
equator because surface conduction is driven by the tan- 
gential component of the bulk electric field (see Figure|SJ). 
Since the double layer in the northern hemisphere is pre- 
dominantly anions, the surface ion transport is from the 
poles towards the equator. A similar argument in the 
southern hemisphere shows that surface transport of the 
majority ion is again from the poles towards the equator. 
Thus, influx of ions into the double layer occurs in the 
regions near the poles and outflux occurs by the equator. 

The dominance of a single species within the double 
layer for each hemisphere leads to a small conundrum: 
how does the bulk electrolyte near the surface of the 
sphere remain locally electroneutral? In the northern 
hemisphere, it would seem that more anion should be 
absorbed at the pole and ejected at the equator leading 
to bulk charge imbalance. Analogous reasoning involving 
cation leads to the same conclusion in the southern hemi- 
sphere. The resolution to the conundrum comes from re- 
membering that both diffusion and electromigration con- 
tribute to transport. The imbalance in the normal flux 
required to sustain an imbalanced concentration profile 
in the double layer is achieved by carefully balancing dif- 
fusion (which drives both species in the same direction) 
and electromigration (which drives the two species in op- 
posite directions) so that the normal flux of the appropri- 
ate species dominates. For example, at the north pole, 
an electric field pointing away from the surface of the 
sphere will suppress the cation normal flux towards the 
surface while enhancing the anion normal flux. 

In this situation, the electric field plays a similar role 
as the diffusion potential for electrochemical transport in 
an electroneutral solution. Recall that when the cation 
and anion have different diffusivities, the electric field 
acts to slow down the species with the higher diffusiv- 
ity and speed up the species with the lower diffusivity in 
such a way that both species have equal flux densities. 
In the current situation, the electric field serves to cre- 
ate the necessary imbalance in the cation and anion flux 



densities so that the surface excess concentrations within 
the double layer can be maintained while preserving local 
electroneutrality in the bulk. 



VI. LINEAR RELAXATION FOR ARBITRARY 
DOUBLE LAYER THICKNESS 

A. Debye-Falkenhagen Equation 

Although the focus of this paper is on nonlinear relax- 
ation, leading to the steady state described in the pre- 
vious section, it is instructive to first consider linear re- 
sponse to a weak applied field, where exact solutions are 
possible. Moreover, the linear analysis also has relevance 
for the early stages of nonlinear relaxation before sig- 
nificant double layer charge has accumulated, such that 
max{e(#)} <C kT/e, as long as the metal surface is un- 
charged when the field is applied. The linear response 
also allows a more general analysis, including AC peri- 
odic response, in addition to our model problem with 
sudden DC forcing. 

When the potential drop across the particle is much 
smaller than the thermal voltage (E a <C 1), it is possible 
to analyze the response of the system without assuming 
that the double layers are thin; that is, we need not as- 
sume that e is small and may describe the system using 
the full unsteady PNP equations, {EI) - Instead, 
we assume that the response of the system is only a small 
deviation from the equilibrium solution: 



1 , p = , <j> — Er cos 9. 



(74) 



In this limit, we can linearize the ionic concentra- 
tions around a uniform concentration profile so that 
c = 1 + (5c. Using this expression in (|2D|> and mak- 
ing use of Poisson's equation (|13J) to eliminate the elec- 
tric potential, we find that the charge density, p = 
(c + — c_) /2 = (Sc + — £c_) 1 2. obeys the (dimensionless) 
Debye-Falkenhagen equation |6(| : 



dp 
at 



9 i 



(75) 



Similarly, the flux boundary condition corresponding to 
this equation reduces to 



dn dn 



(76) 



Note that (|75|l and Poisson's equation (together with the 
no-flux and Stern boundary conditions) are a linear sys- 
tem of partial differential equations. Thus, we can take 
advantage of integral transform techniques. 



B. Transform Solutions for Arbitrary e and S 

Since for weak applied fields, the model problem is a 
linear, initial value problem, it is natural to carry out the 
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analysis using Laplace transforms in time. Transforming 
the Debye-Falkenhagen and Poisson equations, we obtain 



V 2 p 
-e 2 V 2 



I3 2 p 
P 



where 



f3(sf 



(77) 
(78) 



(79) 



and the check accents 
transformed variables, 
tions become 



(") are used to denote Laplace 
Similarly, the boundary condi- 



dp_ + 80 = Q 
dn dn 

(p + oe— = vs 
on 

-V0 — > EqS^ 1 as r — > oo. 



(80) 

(81) 
(82) 



To solve the resulting boundary value problem, we take 
advantage of the spherical geometry to write the solu- 
tion in terms of spherical harmonics. Since p satisfies 
the modified Hclmholtz equation, we can expand it in a 
series with terms that are products of spherical harmon- 
ics, Y[ m (6,(f>), and modified spherical Bessel functions, 
ki((3r). Moreover, we can reduce the series to a single 
term 

p{r,6,4>)=Rk 1 {fir) Y? {6,0) = Rid (fir) cos 6 (83) 

by taking into account the symmetries of the charge den- 
sity: (i) axisymmetry, (ii) antisymmetry with respect to 
9 = 7r/2, and (iii) the dipolar nature of the externally 
applied field. Note that we have only retained the term 
involving the modified spherical Bessel functions that de- 
cays as r — > because p vanishes at infinity. Similarly, 
the general solution for <f> may be expressed as 

4>(r,0,<j>) = -E s- 1 rcos6+A+B^- + C k x (fir) cos 9 

.(84) 

where the first term accounts for the boundary condition 
on the electric field at infinity, the next two terms are 
the general solution to Laplace's equation possessing the 
required symmetries, and the last term is the particular 
solution to Poisson's equation. Note that we have left 
out the monopolar term in the potential because it is 
only necessary for charged spheres. Our analysis is not 
made any less general by neglecting this term; the case 
of a charged sphere in an weak applied field is handled 
by treating it as the superposition of a charged sphere in 
the absence of an applied field with an uncharged sphere 
subjected to an applied field. 

The coefficients in (|83fl and 11841 are determined by en- 
forcing Poisson's equation and the boundary conditions 
(|SU|> - (jSTJ). Plugging and (JSIJl into Poisson's equa- 
tion lO, we obtain 



C = 



R 



(85) 



To apply the boundary conditions, note that on the 
— I 

dr I r— 1 



sphere -§-- = — -Hz L , so that 



= EoS^ 1 cos 9 + 2B cos 9 1 ' 



dn 



(fie) 2 dr 



(86) 



where the last term was obtained by using the relation 
between C and R. Thus, the no-flux boundary condition 
(|8"Uj) becomes 







E n s 1 cos( 



E n s- 



2B cos 9 



2B 




dp 

dr 



r=l 



(87) 



Similarly, the Stern boundary condition (|81|) becomes 

- 1 = A + 

(5e - 1) Eos- 1 + (1 + 2Se) B 



vs 



(0* 



By independently equating the coefficients of the differ- 
ent spherical harmonics in (|87|l and (|88|l , we obtain (after 
a little algebra) 



A = vs 
R 



-3EOS- 1 (Pef 



2fci(/3)+ 1 - (fie) 2 (l + 26e) ^[((3) 



(89) 
(90) 



B = — 



E oS - 



R 



-1 )pk' 1 ((3). (91) 



Finally, by writing k\(x) in terms of elementary func- 
tions 



ki(x) 
we can express R as 

R=- 



e- x (x + 1) 



SEnS' 1 



-0 



(i + ! + r>)u + 2^)-7jV 



(92) 



(93) 



Following Bazant, Thornton, and Ajdari 0, we focus 
on times that are long relative to the Debye time (t = 
0(e 2 ) in dimensionless units). In this limit, s -C 1/e 2 so 
that the charge density on the surface of the sphere is 
given by 



p(r = 1,0,8) 



-Kps- 1 

1 + T p S 



cosf 



with 



7 



3£ Q (l + 6) 
27 

(l + 2Se)(l + e)-Se 



2 7 (l + e) 
(1 + 26e) (1 + e) + S. 



(94) 



(95) 

(96) 
(97) 
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Inverting the Laplace transform, we see that at long 
times, the charge at the surface of the sphere has an 
exponential relaxation with a characteristic time on the 
order of e: 



p{r = l,6,t) ~ -K p (l- e-^*) 



cos 9. 



(98) 



Note that r p is on the order of the dimensionless RC time, 
e, which is much larger than e 2 , the dimensionless Debye 
time. 

To obtain the linear response of cylinders, we follow 
the same procedure as above. The main differences are 
that the series solution is written in terms of a cosine se- 
ries with the radial dependence given by modified Bessel 
functions. Without going through the algebra, the re- 
sults for cylinders are given in Table IIIII along side the 
analogous results for spheres. 



-K q (l - e- l / T ") cosO with 
3E e 



K q = 



7 



(l + 2$e)(l + e) 



27 



(101) 
(102) 



As with p, we see that the characteristic relaxation time 
for q is on the order of the dimensionless RC time. 



E. Time Scales for Linear Response 

We have seen that at long times both p and q relax ex- 
ponentially with characteristic time scales on the order 
of the RC time, e. However, as Figure IT51 shows, the re- 
laxation times for the two quantities are not exactly the 
same and have a nontrivial dependence on the diffuse 
layer thickness, e, and Stern capacitance, S. Notice that 



C. Response to a Weak, Oscillatory Field 

Due to the close relationship between Fourier and 
Laplace transforms, the algebra involved in analyzing the 
response of the sphere to a weak, oscillatory field is al- 
most identical to the response to a suddenly applied field. 
Thus, for sufficiently low frequencies (ui <C 1/e 2 ), we can 
immediately write down the response to a weak, oscilla- 
tory field of the form E = E Q Re (e luJt ) : 



p(r = l,9,t) = -K p Re 



1 + iCJTr 



cos 9. (99) 



D. Accumulated Surface Charge Density 

Because of its importance in many physical processes, 
the accumulated surface charge density, q, is an interest- 
ing quantity to consider. It is easily computed from the 
(volume) charge density p by integrating it in the radial 
direction from r = 1 to r — oo. While the identifica- 
tion of this integral with a surface charge density really 
only makes sense in the thin-double layer limit, the ac- 
cumulated surface charge density is still worth examin- 
ing. Fortunately, the integral is straightforward because 
k±(z) = — fcg(z) [6l| and the radial dependence of the 
charge density is independent of the angle. Using these 
observations, the surface charge density is given by 



R k Q (P) 

q(Q,s) = — -g — cos( 

Re-? 
= — cos*. 



(100) 



In the long time limit s <C 1/e 2 , we find that the sur- 
face charge density an exponential relaxation q{9, t) = 




FIG. 13: Exponential relaxation time constants for 
the charge density p and the accumulated surface 
charge density q at weak applied fields as a function 
of e and 8. The left panel shows the relaxation time 
constant for the charge density at the surface of the 
sphere, r p (r = 1). The right panel shows the relax- 
ation time constant for the accumulated surface charge 
density, r q . 



for infinitely thin double layers (e = 0) the relaxation 
times for the surface charge density and the accumulated 
charge density are identical. This behavior is expected 
since for thin double layers, almost all of the charge den- 
sity in the diffuse layer is located very close to the surface 
of the particle. In this limit, the relaxation time has a 
strong dependence on the Stern capacitance. For thick 
double layers (e> 1), this dependence on the Stern ca- 
pacitance disappears and the relaxation time curves for 
all 5 values converge. 

Physically, the difference in the relaxation times for the 
surface charge and the accumulated charge densities for 
nonzero e values is an indication of the complex spatio- 
temporal structure of the double layer charging. For thin 
double layers, the difference r p and r q is relatively small 
because the charge in the double layer is restricted to a 
thin region. For thick double layers, however, the dif- 
ference in relaxation times is accentuated because the 
charge in the double layer is spread out over a larger spa- 
tial region, which does not necessarily charge uniformly. 
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TABLE III: Table of formulae for the linear response of metallic cylinders and spheres to weak applied electric fields. 



Cylinder 



Sphere 



R 
B 

Kn 



K 



RKx{j3r) cos (9 
- 1 -£ oS - 1 rcose+ 



03 <0 2 



-2E s~ 1 (/3e) 2 



-E oS - x -R^-l)f3Km 

2E a K 1 (l/t) 
K^l/^-SK'^l/e) 

2[ifi(l/ £ )-Xj(l/ E )] 
2E o eK (l/e) 



Rk!{/3r) cos 6 
vs- 1 - E oS - l r cos ^iLM] c< 

3Eo(l+6) 
2[(l+2«e)(l+e)+c5] 

f (l+2^)(l+ £ )-j £ 1 
fc \ 2(l+e)[(l+2«e)(l+e)+,5] J 

3B„t 



ifl(l/e)-«X((l/ E ) 



( e +^Tw) Kl(1/£)+ [ 1+2oe - 5 5§Tw]^ (1/e)+,5x " (1/E) 

2[K 1 (l/ f )-Ki(lA)J 



(l+2ie)(l+e)+« 



. f (1+236) (1+e) 1 
\ 2[(l+24e)(l+e)+«] J 



a measured from vertical axis. 



In fact, that r p (r = 1) is smaller than r g for larger values 
of e fFigurc lT3|) suggests that regions closer to the surface 
of the sphere charge faster than regions that are further 
away. 



VII. WEAKLY NONLINEAR RELAXATION 
FOR THIN DOUBLE LAYERS 

A. Dynamical Regimes in Space-Time 

For weak applied fields in linear response, the compli- 
cated dependence of the Laplace transform solution for 
large s (short times) above hints at the presence of multi- 
ple time scales in the charging dynamics. In this section, 
using boundary-layer methods, we derive asymptotic so- 
lutions in the thin double-layer limit (e <C 1) at some- 
what larger electric fields (defined below) for the con- 
centrations and electric potential by solving the leading 
order equations at the two dominant time scales: (1) the 
RC time Xna/D and (2) the bulk diffusion time a 2 /D. 
We proceed by seeking the leading order term (and in 
some cases, the first-order correction) to the governing 
equations (|19fl - H21|) with both the spatial and the time 
coordinate scaled to focus on the space-time region of 
interest. 

For the analysis, it is important to realize that the 
space-time domain is divided into five asymptotically dis- 
tinct regions (see Figure H3) l. At the RC time, there ex- 
ist three spatially significant regions: (i) an 0(e) quasi- 
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FIG. 14: Five asymptotically distinct regions of space- 
time that govern the dynamic response of a metal col- 
loid sphere to an applied electric field. Note the nested 
spatial boundary layers at the RC time (t = O(e)). 



steady double layer, (ii) an 0(y/e) dynamically active 
diffusion layer, and (iii) a quasi-equilibrium, uniform 
bulk electrolyte layer with time-varying harmonic elec- 
tric potential. At this time scale, the charging dynamics 
are completely driven by the 0(y/e) diffusion layer. At 
the diffusion time, there are only two important spatial 
regimes: (i) a quasi-steady double layer and (ii) a dy- 
namic bulk that evolves through locally electroneutral, 
diffusion processes. 
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B. Dynamics at the RC Time 

1. Uniform Bulk and Equilibrium Double Layers 

To examine the dynamics at the RC time, we rewrite 
(|19[) - Q21[l by rescaling time to the RC time using t = t/e: 



8c 

dt 
dp 

dt 



eV ■ (Vc ■ 



eV ■ (Vp- 



-cV0) 



-e 2 V 2 (/. = p. 



(103) 

(104) 
(105) 



Since the spatial coordinate is scaled to the bulk length, 
the leading order solution of these equations describes the 
dynamics of the bulk during the double layer charging 
phase. Substituting a regular asymptotic expansions of 
the form 



c(x, t) ~ c + eci + e 2 c 2 + 



(106) 



into equations (|103fl - l|105fl yields a hierarchy of partial 
differential equations. By sequentially solving the equa- 
tions using the initial conditions, it is easy to show that 
the "outer" solutions at the RC charging time scale are 



c (x,t) 

Cj(x,t) 

Pj(x,t) 



= 1 



0, j — 1, 2, 3, 
0, j — 0, 1, 2, 



(107) 
(108) 
(109) 



with is harmonic at all orders. In other words, the bulk 
solution can be completely expressed ( without a series ex- 
pansion) as a uniform concentration profile, c (x,t) = 1, 
and a time-varying harmonic electric potential, <p. Note 
that by taking advantage of spherical geometry and ax- 
isymmetry in our problem, we can write the potential as 
a series in spherical harmonics with zero zonal wavenum- 
ber (i.e., Legendre polynomials in cos6>): 



r,9,t) = -E Q r cos 9 + ^2 Pi (cos 9) 



Ai(t) 
r i+i 



(110) 



where the radial dependence of each term has been se- 
lected so that 4> automatically satisfies Laplace's equation 
at all times. 

At the RC charging time, the 0(e) double layer is in 
quasi-equilibrium, which is easily verified by rescaling the 
spatial coordinate normal to the particle surface by the 
dimensionless Debye length e. As mentioned earlier, a 
quasi-equilibrium double layer possesses a structure de- 
scribed by the Gouy-Chapman-Stern model. For conve- 
nience, we repeat the GCS solution here: 



c± = 



c = ccosh-0 , p = — csmhijj 



i>(z) 



4 tanh 



1 (tanh(C/4)e 



— V&2 



(111) 



Recall that ip is the excess voltage relative to the bulk, 
tp = 4> — cf> and that the tilde and hat accents are used to 



indicate quantities within and just outside of the double 
layer, respectively. Also, we reiterate that unlike the ID 
case, £ is not a constant but is a function of spatial along 
the electrode surface. 

Since the bulk concentration at the RC time is uni- 
form, the double layer structure is given by l|lllfl with c 
set equal to 1. Notice that at the RC time, bulk concen- 
tration gradients have not yet had time to form, so the 
variation of the diffuse layer concentration and charge 
density along the surface of the sphere is solely due to a 
non-uniform zeta-potential. 



2. O (\/e) Diffusion Layer 

The analysis in the previous section leads us to an ap- 
parent paradox. The dynamics of the system seem to 
have been lost since both the bulk and the boundary 
layers are in quasi-equilibrium at leading order; neither 
layer drives its own dynamics. As discussed in [l|, the 
resolution to this paradox lies in the time-dependent flux 
matching between the bulk and the boundary layer. Un- 
fortunately, it is inconsistent to directly match the bulk 
to the boundary layer; there must exist a nested O (i/e) 
diffusion layer in order to account for the build up of both 
surface charge and surface excess neutral salt concentra- 
tion. 

Mathematically, the presence of the diffusion layer at 
the RC time scale appears as a dominant balance in the 
transport equations by rescaling the spatial coordinate 
in the normal direction to the surface by ^ft to obtain: 



dc 
dt 



e 2 V s • (V s c + pV 5 



+ 



d f dc _d(f> 



dz' \dz' ^ dz' 



^ = e 2 V s • (V s p + cV s 
dt K F 



dz' V dz' dz' 



-e 2 V 2 



dz'' 



(112) 

(113) 
(114) 



Here the bar accent denotes the "diffusion layer" solution 
at the RC time and z' = Z/ ^/e is the spatial coordinate in 
the direction normal to the surface. Notice that at this 
length scale, the system is not in quasi-equilibrium as 
it is at the bulk and Debye length scales. It is, however, 
locally electroneutral at leading order as a result of (|114fl . 

As the double layer charges, it absorbs an 0(e) amount 
of charge and neutral salt from the 0(y/e) diffusion layer. 
Therefore, we expect that concentration changes within 
the diffusion to be on the order of which motivates 
the use of an asymptotic expansion of the form 



c(x, t) ~ c + 



£ 1/2 C 1/2 



e d 



(115) 
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Using this expansion in 1|112|) 
leading order equations are: 



d_ 
~dz' 



cq- 



dc 
dt 

00 



pi4p . we find that the 
d 2 c 



dz' 2 
= 



Po = 0. 



(116) 

(117) 
(118) 



The initial conditions and boundary condition as z' — > oo 
for these equations are Go(t = Q) = 1 and cq(z' — > oo) = 
1, respectively. The boundary condition at z' = is given 
by flux matching with the double layer. Rescaling space 
and time in (|60|1 and (|61|l . we find that the appropriate 
flux boundary conditions for the diffusion layer are 



dt 

dw 



eV s • (qV s logc + wV s 



dt 



Y = eV s • (fflV s logc + qV s 



1 
1 



dej) 



dz 

dc 



7 (H9) 



d: 



-. (120) 



Thus, the leading order flux boundary conditions, which 
appear at 0(l/y/e), are 



Co 



cm 

dz' 
dc Q 
dz' 



(121) 
(122) 



The leading order solutions in the diffusion layer, Cq = 1 
and <fio = 4> (Z — > 0), are easily determined by applying 
the initial and boundary conditions to 1)116(1 - 1118(1 . 

To obtain dynamics, we must examine the first-order 
correction to the solution. At the next higher order, the 
governing equations are 



dCj/2 

dt 

h/2 



"J Cl/2 

dz' 2 



d 2 



dz' 2 

Pl/2 



(123) 

(124) 
(125) 



where we have made use of the leading order solution to 
simplify ((124(1 . Again, the initial conditions and bound- 
ary condition at infinity are simple: c-y^it — 0) = 
and Ci/2{z' oo) = 0. The flux boundary conditions 
at z' = 0, however, are more interesting because they 
involve the charging of the double layer: 



(126) 
(127) 



dq 


O01/2 


dt 


dz' 


dw 


dci/2 


~dt 


dz' 



2 = 



Thus, simple diffusion of neutral salt at x = 0(i/e) driven 
by absorption into the 0(e) double layer dominates the 
dynamics of the diffusion layer. From ((124(1 and ((126(1 . 
we see that the electric potential possesses a linear profile 



with slope given by rate of surface charge build up in the 
double layer: 



(Z^0) + e 



1/2 



dt 



(128) 



Note that constant term at 0(y/e) is forced to be zero 
by matching with the bulk electric potential since all 
higher order corrections to the bulk potential are identi- 
cally zero. 



3. Effective Boundary Conditions Across Entire Diffusion 
Layer- 
It is precisely the fact that the electric potential has 
the form 1(128(1 that justifies the common approach of 
asymptotic matching directly between the bulk and the 
double layer [ll[iall3- For instance, the time-dependent 
matching used by Bazant, Thornton, and Ajdari [l| rests 
on the implicit assumption that the leading order electric 
field in the diffusion layer is a constant and appears at 
0{\/~i) so that 



dZ 



1 



z=o 



d<f>_ 
~dz~' 



?l/2 



dz' 



^oo dt 

(129) 

Similarly, the definition of the zeta-potential and the 
Stern boundary condition ((31) I) across the entire diffusion 
layer require that <p ~ <p (Z — >0)at leading order. 



4- Leading-order Dynamics 

Using the results of our discussion, we now derive the 
leading order equations that govern the charging dynam- 
ics on the surface of the sphere. Since charging is non- 
uniform over the electrode surface, the equations take 
the form of partial differential equations on the surface 
of the sphere. Defining \E f (6 l , 0) = v — (j>(r — 1,0, (f>), 
we can write the Stern boundary condition ((30(1 and the 
double layer charging equation ((129(1 as 



C + 2<5sinh(C/2) = f 



9* 



dn 



(130) 
(131) 



where we have introduced the leading order differential 
double layer capacitance C(4') = —dq/d^ and used the 
fact that c = 1 at the RC time. Together with (J2U, fT3nj) 
and ((131(1 form a complete set of boundary conditions for 
the leading order electric potential in the bulk region. To 
compute the double layer capacitance, we can combine 
(ISUll with (|HU|> to obtain 



C = 



1 



sech(C/2) 



(132) 



Since C depends on ^ (via £), the charging equation 
((131(1 is nonlinear which makes the problem analytically 
intractable. 
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For small applied fields (where £ « E Q is a reasonable 
approximation), analytical progress can be made by lin- 
earizing the the double-layer capacitance around ( = 
to obtain C ~ 1/(1 + 5). This approximation leads to 
a linear charging equation making it possible to solve 
the equations as an expansion in spherical harmonics. 
Substituting the expansion (| 1 1 0|) into i|131[l and the def- 
inition of '3/ , we obtain a decoupled system of ordinary 
differential equations for the expansion coefficients: 



dAo 
dt 



^ + 2(1 + S)A 1 
at 

A A. 

-l + (l + l)(l + 6)A l 



dv 
dt 



= , l> 1. 



(133) 



To retain generality, we have allowed for the possibility 
that the applied electric field and surface potential are 
time-varying. For the case of a steady surface potential 
and uniform applied field, we find that the bulk electric 
potential is 



V< 1+a > f -£„rcos< 



1 



2r 3 ( 



1 -3e 



-2(l+5)t 



(134) 

where we have assumed that both the surface potential 
and the electric field are suddenly switched on at t = 
0. The initial condition in this situation is that of a 
conducting sphere at potential v in a uniform applied 
electric field E Q : (j) = v — E o rcos0 (l — l/?" 3 )- 

Using (|134|) , the double layer potential and total diffuse 
charge are easily determined to be 



$ = 



■( 



-(l+<5)* 



£ o cos0(l-e- 2 ( 1+ 1 



') 



r 



1 - e- [1 +^ 



2(1 + 5) 



E D cos ( 1 — e 



-2(l+<5)t 



(135) 



(136) 



These results are consistent with the calculations by 
Bazant and Squires [itiL Il7j , which is expected since the 



low applied field limit corresponds to the regime where 
the total diffuse charge is linearly related to the zeta- 
potential (and therefore the total double layer poten- 
tial). It is worth mentioning that in the common situa- 
tion where the surface potential is set well before t = 0, 
then the first term in each of the above three expressions 
is not present. Analogous double layer charging formulae 
for the cylinders are shown in Table llVl 



5. Numerical Model 

For larger applied fields, the nonlinear double layer 
capacitance forces us to use numerical solutions to gain 
physical insight. Since the bulk electric potential is a 
time- varying harmonic function, it is most natural to nu- 
merically model the evolution equation for the potential 
using a multipole expansion with harmonic terms. Trun- 
cating Qll()|l after a finite number of terms yields a dis- 
crete solution of the form: 



ct>(r,6j) 



-E n rcost 



' J y 

1=0 



^P(cos0) (137) 



where the unknowns are the time-dependent coefficients 
in the expansion. By using (|137fl and enforcing that l|131|) 
is satisfied at the collocation points B L = in/N, we de- 
rive a system of ordinary differential equations for the 
coefficients Ai(t): 



CP d 4 

dt 



E P(:,2) + QA 



(138) 



where A is the vector of expansion coefficients, P and Q 
are collocation matrices that relate the expansion coeffi- 
cients to <fi and d(j)/dn, respectively, and C is a diagonal 
matrix that represents the double layer capacitance at 
the collocation points. Note that P(:,2) (which repre- 
sents the second column of P in MATLAB notation) is 
the discrete representation of Pi (cos 6) — cos 9, so the 
term £' P(:,2) accounts for the applied background po- 
tential. The explicit forms for the collocation matrices 
and the discrete double layer capacitance are given by 



P (cos#o) Pi (cos 9 Q ) 
P (cos6>i) Pi (cos 0i ) 



Pv(cos O ) 
P/v(cos0i) 



Po(cos0at) Pi(cos0at) ... Pn{cos9n) 



Q 



P o (cos0 o ) 2Pi(cos0 o ) 
P o (cos0i) 2Pi(cos0i) 



(A+l)Pv(cos0 o )' 
(A+l)Pv(cos0i) 



P o (cos0at) 2Pi(cos0at) ... {N + l)P N (cos9 N ) 
C = diag(C7(*(0o) ) C(*(0i),...,C(*(0 J v)). (139) 



The system of ODEs for the expansion coefficients is eas- tiplying (|138|l through by (CP) 1 and writing a simple 
ily solved using MATLAB 's built-in ODE solvers by mul- 
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TABLE IV: Table of formulae for double layer charging of metallic cylinders and spheres at weak applied electric fields at the 
RC-time. In these formulae, the potential, v, of the particle is set to zero. 



Cylinder 11 



Sphere 



-£ o rcos0[l + 4 r (l-2e-( 1 +^ t )] 



q 



2E cos 6 1 - e 



-(l+S)t 



-jj^Eo cos6U 1 — e 



-(l+5)t 



-£ rcose[l + ^(l-3e- 2(1+,s)t 
§£ cos6i (l - e" 2( 



-2(l+«)t 



a 9 measured from vertical axis. 




v = 
E = 5 
8 = 0.1 



4 ~~ 2 4 6~ 8 10 
t 

FIG. 15: Time evolution of the dominant coefficients 
in the Legendre polynomial expansion of the bulk elec- 
tric potential in the weakly nonlinear regime at the RC 
time. In this figure v = 0, E — 5 and 5 — 0.1. Notice 
that the dipolar term dominates the solution. 




4 ~~ 2 4 6~ 8 10 
t 



FIG. 16: Time evolution of the dominant coefficients 
in the Legendre polynomial expansion of the bulk elec- 
tric potential in the weakly nonlinear regime at the RC 
time when the sphere has a nonzero applied voltage. 
In this figure v = 3, E — 5 and 8 = 0.1. Note that 
both symmetric and antisymmetric terms make non- 
negligible contributions. 



function to evaluate the resulting right-hand side func- 
tion. 



6. Dipole-Dominated Charging 

From the numerical solution of the charging equa- 
tion (|131|1 . we see that when the sphere is electrically 
grounded, charging is dominated by the dipolar contribu- 
tion to the response (see Figure lT5)l . While the nonlinear 
capacitance does in fact allow higher harmonics to con- 
tribute to the response, the higher harmonics only play 
a small role even at larger applied fields. As expected, 
when the sphere is kept at zero voltage, the antisymme- 
try between the upper and lower hemisphere is not bro- 
ken and only odd terms contribute to the series solution 
(|110fl . However, as shown in Figure [TBI if a nonzero po- 
tential is applied to the sphere, all harmonics contribute 
to the solution. In this case, the dominant contributions 
come from the monopolar and dipolar terms. 

The dipolar nature of double layer charging forms the 
foundation of much of the work on the charging of col- 
loid particles over the past half century. For instance, 
the non-equilibrium double layer is often characterized 



in terms of the induced dipole moment |56j. As shown 
in (I134[) - (|136[) . the monopolc and dipole contributions 
are the only contributions in a linearized theory. Our 
numerical investigations demonstrate that, even for the 
nonlinear theory, the monopolc and dipole response dom- 
inates the total response. Our results provides additional 
theoretical support for the focus on the dipole response 
when studying colloid particles in applied fields. 

7. Extended Double Layer Charging 

The slowing of double layer charging is one important 
consequence of nonlinearity in the double layer capaci- 
tance (see Figure IT7f) . However, as shown in Figure IT7I 
extended charging only occurs for 8 <C 1. Mathemati- 
cally, we only see slowed charging at small values of 5 be- 
cause the double layer capacitance 1)1321) can only become 
as small as 1/6, which occurs at large zeta-potentials. 
For sufficiently small 6, charging is slowed at higher ap- 
plied fields because the sech(£/2) term in the denomina- 
tor of the double layer capacitance (|132f) becomes negli- 
gible when C > -21n(<5/2). 
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v = 
E = 5 
5 = 0.01 



1 2 3 4 5 5 10 15 20 25 

t t 

FIG. 17: Time evolution of the dominant coefficients 
in the Legendre polynomial expansion of the bulk elec- 
tric potential in the weakly nonlinear regime at the 
RC time for low (left) and high (right) Stern capaci- 
tance values. In these figures v = and E — 5. Note 
that double layer charging is retarded when 5 is small 
but that the slowdown in double layer charging is sup- 
pressed for larger 8 values. 



Notice that there is no time dependence for either the 
concentration or the electric potential. It is worth men- 
tioning that for a general geometry, the initial condition 
is a uniform concentration profile with the electric poten- 
tial of an insulator in an applied field and the boundary 
conditions (which are consistent with the initial condi- 
tions) are 

c —- = and -5— = 0. (146) 
On on 

At the diffusion time scale, the double layer contin- 
ues to remain in quasi-equilibrium, so its leading order 
structure is given by 1)111(1 with the bulk concentration 
set equal to 1. However, unlike the double layer at the RC 
time scale, the leading order zeta-potential is not evolv- 
ing, so the leading order double layer structure is static 
in time. 



C. Dynamics at the Diffusion Timescale 



2. Higher-Order Bulk Diffusion 



In this section, we examine the response of the system 
at the diffusion time scale. We find that the only dy- 
namic process is diffusion of neutral salt within the bulk 
in response to surface adsorption that occurred at the 
RC time scale. Since the total amount of neutral salt ab- 
sorbed by the diffuse charge layers during the charging 
phase is an O(e) quantity, the bulk concentration only 
needs to decrease by O(e) to compensate. Thus, we find 
that dynamics are not present at leading order; rather, 
they appear only in higher-order corrections. Also, at 
the diffusion time, surface transport, which is negligible 
at the RC time scale, becomes important. 



1. Leading Order Bulk and Double Layer Solutions 

Substituting an asymptotic series into (|19J) - i|21|l . we 
obtain the leading order bulk equations: 



(140) 



V 2 (/> = 0. (141) 

Applying the initial conditions (obtained by matching to 
the solution at the RC time) and the leading order bound- 
ary conditions (derived from the effective flux boundary 
conditions (|60|l and Hfilfl ) yields the simple leading order 
solutions for a sphere 



Co(x,t) = 1 

(j>o(x,t) = -E Q r cos 6 ^1 + ^3 

and a cylinder 

c (x,t) = 1 

(f>o(x,t) = -E o rcos0 (l + 



(142) 
(143) 

(144) 
(145) 



In order to see dynamics, we need to consider the first- 
order correction to (|140fl and l|141fl : 



dc i 



V 2 ci 



V 2 0i = -Vci • V0o- 



(147) 
(148) 



The boundary conditions for these equations are a bit 
more complicated. Using ijrjU)) - ijfTTl) and taking into ac- 
count the leading order solutions, we find that the bound- 
ary conditions for the 0(e) equations are 

q 6 + (t) = V s • (woV^q) - ^p- (149) 

cm 

w S+(t) = V s • (goV^o) - (150) 

on 

where go an d wq are the leading order equilibrium sur- 
face charge density and surface excess neutral salt con- 
centration. As mentioned earlier, at the diffusion time 
scale, these quantities are static in time. Also, note the 
presence of the delta-functions in time, which account 
for the "instantaneous" adsorption of charge and neutral 
salt from the bulk during the at the RC-time |l|. 

Mathematically, the appearance of the delta-functions 
is a consequence of the connection between the time 
derivative of double layer quantities at the two time scales 
in the asymptotic limit e — > 0. To illustrate this connec- 
tion, consider the time derivative of the surface charge 
density, q. Let t and t be scaled to the diffusion and RC 
times, respectively, so that t = et. At these two time 
scales, the surface charge density can be written as q(t) 
and q (t) which are simply related by q(t) — q (t) . The 
time derivatives, however, are related by 



dt 



(*) 



1 dq 



e dt 



Ht/e) 



(151) 
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Therefore, for nonzero t, ^| = in the asymptotic limit 
because || (t/e) approaches zero faster than linearly as 
e — > (as an example, see (I136[) ). In contrast, for t = 0, 
^| is infinite because ^f(O) has a fixed nonzero value. 
Next, consider the following integral with ti > 0: 

It 3* = i A = * " $ ^ ' (152) 

In the asymptotic limit, the integral approaches zero for 
nonzero t\ but approaches <z(oo) when t\ equals zero. 
Putting the above properties together, we see that ^| (i) 
is indeed a one-sided delta-function of strength q(co) — 

In contrast to one-dimensional systems, the boundary 
layers play a more active role in the evolution of the bulk 
concentrations because surface conduction continues to 
play a role well beyond the initial injection of ions at 
t = 0. Note, however, that the surface conduction terms 
in 1149|) and l|15U|) are static, so they essentially impose 
fixed normal flux boundary conditions on the 0(e) bulk 
equations. 



3. Comparison with the Analysis of Dukhin and Shilov 

It is interesting to compare and contrast our weakly 
nonlinear analysis with the work of Dukhin and Shilov 
[42l |45| on the polarization of the double layer for highly 
charged, spherical particles in weak applied fields. In 
both cases, bulk concentration variations and diffusion 
processes appear as a higher-order correction to a uni- 
form background concentration and are driven by surface 
conduction. However, the significance of the surface con- 
duction term arises for different reasons. As mentioned 
earlier, the small parameters that controls the size of the 
correction is different in the two analyses. In Dukhin 
and Shilov's analysis, the small parameters are e and E a . 
Because they essentially use asymptotic series in E D as 
the basis for their analysis, they require that the double 
layer be highly charged in order for surface conduction to 
be of the same order of magnitude as the 0{E ) normal 
flux of ions from the bulk. In other words, because the 
size of the surface conduction is proportional to E a , in 
order for surface conduction to be of the same order of 
magnitude as the normal flux, the surface charge density 
must be an 0(1) quantity: eq = 2esinh(Co/2) = 0(1). 
In contrast, we use asymptotic series in e in our analysis 
and do not restrict E Q to be small, so the surface conduc- 
tion and normal flux of ions from the bulk have the same 
order of magnitude for small O(e) surface charge densi- 
ties regardless of the strength of the applied electric field 
(as long as it is not so large that the asymptotic anal- 
ysis breaks down). Thus, our weakly nonlinear analysis 
complements the work of Dukhin and Shilov by extend- 
ing their analysis to stronger applied electric fields and 
the case where the surface charge density is induced by 



the applied field rather than fixed by surface chemistry 
of the colloid particle. 



VIII. STRONGLY NONLINEAR RELAXATION 

A. Definition of "Strongly Nonlinear" 

The weakly nonlinear analysis of the thin double-layer 
limit in the previous section assumes that the dimension- 
less surface charge density a and surface salt density j3 
remain uniformly small. For our PNP model, we can 
estimate when this assumption becomes significantly vi- 
olated using Eq. with Q = E a/ (1 + S) and C = C , 
which occurs at fields above a critical value at least a few 
times the thermal voltage, 

*o>^(l + ^W(f). (153) 

As discussed in section IIVI this condition also implies 
that surface adsorption of ions from the bulk is larger 
enough to trigger significant bulk diffusion and surface 
transport through the double layer, in steady steady 
state. However, as noted by Bazant et al. |l|, weakly 
nonlinear asymptotics breaks down during relaxation dy- 
namics at somewhat smaller voltages, a/s/ne < 1, or 
(with units restored) 




(154) 



since large surface adsorption can occur only temporarily 
at certain positions. (The factor 7r in this formula comes 
from a one-dimensional analysis of bulk diffusive relax- 
ation, which does not strictly apply here, but the rough 
scale should be correct. 

Beyond the weakly nonlinear regime, there are two 
main effects that occur: (1) transient, local depletion of 
the leading order bulk concentration and (2) surface con- 
duction at the leading order. As in the case of a steady 
applied field, perhaps the greatest impact of a — 0(1) is 
that we must pay attention to factors of the form ee^ or 
esinh(C), in addition to factors of e, when ordering terms 
in asymptotic expansions. 

In the thin double- layer limit, the boundary layers are 
still in quasi-equilibrium, as long as the bulk concen- 
tration does not approach zero, which would typically 
require Faradaic reactions consuming ions at the con- 
ducting surface [5lL l62l | . Since we only consider ideally 
polarizable conductors, bulk depletion is driven solely by 
adsorption of ions in the diffuse layer, which is unlikely to 
exceed diffusion limitation and produce non-equilibrium 
space charge, although this possibility has been noted pj . 
Therefore, we would like to proceed as in the previous 
sections and treat the bulk as locally electroneutral with 
effective boundary conditions. 
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B. Leading-Order Equations 

Unfortunately, the analysis of the leading order equa- 
tions derived in this manner does not appear to be as 
straightforward as the analysis in the weakly nonlinear 
limit. The main problem is that it seems difficult to de- 
rive a uniformly valid leading order effective boundary 
conditions along the entire surface of the sphere. For 
this reason, we merely present the apparent leading or- 
der equations for the strongly nonlinear regime and leave 
a thorough analysis for future work. 

At the leading order in the bulk, we find the usual 
equations for a neutral binary electrolyte: 

^ = V 2 c (155) 
at 

= V ■ (c o V0 o ) (156) 

with p = 0(e 2 ). The structure of the boundary layers 
is described by GCS theory with the concentration and 
charge density profiles given by 1|> . Effective bound- 
ary conditions for (|155fl _ (|156fl are derived in the same 
manner as for a steady applied field except that unsteady 
terms are retained. Recalling that q and w grow expo- 
nentially with the zeta-potential, we find that both the 
surface conduction terms and the time derivatives of the 
total diffuse charge and excess concentration appear in 
the leading order boundary conditions: 

e-7£- = £ V S • (w \7 s (j) ) - c — (157) 

= eV s • (ftV.fc) - (158) 
at on 

C. Challenges with Strongly Nonlinear Analysis 

It is important to realize that these equations are 
mathematically much more complicated than the anal- 
ogous equations in the weakly nonlinear regime. First, 
the nonlinearity due to the clectromigration term explic- 
itly appears in the bulk equations at leading order; the 
nonlinearity is not removed by the asymptotic analysis. 
Furthermore, the diffuse layer concentrations depend on 
time explicitly through the bulk concentration at the sur- 
face in addition to the zeta potential: 

c±(t) =c±(i)e^W. (159) 

Already, these features of the equations greatly increases 
the challenge in analyzing the response of the system. 

However, the greatest complication to the mathemat- 
ical model in the strongly nonlinear regime is that effec- 
tive boundary conditions l|157fl and (|158fl are not uni- 
formly valid over the surface of the sphere (or cylinder). 
Near the poles, the double layer charges quickly, so the 
time-dependent and surface transport terms in the ef- 
fective boundary conditions become 0(1) at very short 
times. In contrast, the amount of surface charge in the 



double layer near the equator is always a small quantity. 
Thus, it would seem that near the equator, the only sig- 
nificant terms in the effective boundary conditions are 
the normal flux terms. Together, these observations sug- 
gest that the appropriate set of boundary conditions to 
impose on the governing equations l|155|) and (|156fl de- 
pends on the position on the surface of the sphere. More- 
over, the position where the boundary conditions switch 
from one set to the other depends on time as double-layer 
charging progresses from the pole towards the equator. 



IX. CONCLUSION 
A. Predictions of the PNP Model 

In this paper, we have analyzed electrochemical relax- 
ation around ideally polarizable conducting spheres and 
cylinders in response to suddenly applied background 
electric fields, using the standard mathematical model 
of the Poisson-Nernst-Planck (PNP) equations. Unlike 
most (if not all) prior theoretical studies, we have focused 
on the nonlinear response to "large" electric fields, which 
transfer more than a thermal voltage to the double layer 
after charging. We have effectively extended the recent 
nonlinear analysis of Bazant, Thornton and Ajdari [If for 
the one-dimensional charging of parallel-plate blocking 
electrodes, to some new situations in higher dimensions, 
where the potential of the conductor is not controlled. 
Instead, electrochemical relaxation in our model prob- 
lems is driven by a time-dependent background electric 
field applied around the conductor, whose charges are 
completely free to relax. We are not aware of any prior 
analysis of nonlinear response in such problems (whether 
or not using the PNP equations), and yet it is important 
in many applications, in microfluidics, colloids, and elec- 
trochemical systems, where applied fields or voltages are 
often well beyond the linear regime. 

We have focused on the structure and dynamics of the 
double layer and the development of bulk concentration 
gradients. A major goal has been to move beyond the tra- 
ditional circuit models commonly used to study the linear 
response of electrochemical systems. Through a combi- 
nation of analytical and numerical results, we have shown 
that significantly enhanced ion concentration within the 
double layer is a generic feature of nonlinear electrochem- 
ical relaxation around a conductor. By interacting with 
the bulk electric field, the enhanced concentration within 
the double layer leads to large surface current densities. 
In addition, because the double layer does not charge 
uniformly over the surface, tangential concentration gra- 
dients within the double layer itself lead to surface diffu- 
sion. Due to their coupling to bulk transport processes 
via normal fluxes of ions into the double layer, these sur- 
face transport process drive the formation of bulk con- 
centration gradients. 

We have also found that bulk concentration gradients 
are always present to some degree and play an important 
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role in allowing the system to relax to the steady-state. 
For weak applied fields, they are often neglected because 
they only appear as a first-order correction to a uniform 
background concentration profile. However, for strong 
applied fields, the variations in the bulk concentration 
becomes as large as the background concentration, so 
they cannot be ignored. These bulk concentration gra- 
dients lead to bulk diffusion currents which result in net 
circulation of neutral salt in the region near the metal 
colloid particle. 

Another key contribution of this work is a careful 
mathematical derivation of general effective boundary 
conditions for the bulk transport equations in the thin 
double layer limit by applying matched asymptotic ex- 
pansions to the PNP equations. We derive a set of bound- 
ary equations that relate the time evolution of excess sur- 
face (double-layer) concentrations to surface transport 
processes and normal flux from the bulk. An interesting 
feature of the effective boundary conditions is that they 
explicitly involve the small parameter e and may have 
different leading order forms depending on the charac- 
teristic time scale and the magnitude of double layer ion 
concentrations . 

It is beyond the scope of this article, but worth pur- 
suing, to study "strongly nonlinear" dynamics in detail, 
where (according to the PNP equations) the double layer 
adsorbs enough ions to significantly perturb the bulk con- 
centration and cause surface transport to bulk transport 
of ions. The challenging issues discussed at the end of 
Section IVIIII need to be addressed in order to gain a 
deeper understanding of the rich behavior of metal col- 
loid systems in this practically relevant regime. Also, it 
would be beneficial to validate solutions of the effective 
bulk and boundary conditions in the thin double layer 
limit against solutions for the full PNP equations. The 
utility of the thin double layer approximation cannot be 
fully appreciated until this comparison is completed. 



B. The Need for Better Continuum Models 

Unfortunately (or fortunately, depending on one's per- 
spective) , the theoretical challenge of strongly nonlinear 
electrochemical relaxation is much deeper than just solv- 
ing the PNP equations in a difficult regime: It is clear 
that the dynamical equations themselves must be modi- 
fied to better describe the condensation of ions in highly 
charged double layers. As emphasized by Newman |38| . 
the PNP equations are only justified for dilute solu- 
tions, since each ion moves in response to an independent 
stochastic force with constant diffusivity and mobility, 
interacting with other ions only through the mean long- 
range Coulomb force. Important effects in concentrated 
solutions, such as short-range forces, many-body correla- 
tions, steric constraints, solvent chemistry, and nonlinear 
mobility, diffusivity, or permittivity, are all neglected. 

It is tempting to trust the PNP equations in the case of 
a dilute bulk electrolyte around a initially uncharged sur- 



face. However, as we have shown, the model predicts its 
own demise when a large electric field is applied, due to 
enormous increases in surface concentration in the diffuse 
part of the double layer, even if the bulk concentration 
is small. Note that the condition for strongly nonlinear 
relaxation (|I53[) or (|I54|I is similar to the condition for 
the breakdown of the PNP equations - both require ap- 
plied voltages across the double layer only several times 
the thermal voltage. In particular, the Gouy-Chapman 
solution to the PNP equations predicts an absurd con- 
centration of counterions of one per A 3 at the surface for 
a surprisingly small zeta potential around bkT / z + e ~ 0.2 
V even in a fairly dilute 0.01 M electrolyte. This critical 
voltage is routinely applied to the double layer in elec- 
trochemical systems. 

Of course, electrochemists are well aware of this prob- 
lem; in fact, Stern originally proposed the compact layer 
of adsorbed ions, outside the continuum region, as a way 
to cutoff the unbounded double-layer capacitance |63| . 
Since then, various emp irical models of the compact layer 
have appeared [3^, |3||, but steric effects (or other non- 
linearities) in concentrated solutions have received much 
less attention. Borukhov, Andelman and Orland recently 
postulated a continuum free energy for ions, taking into 
account steric repulsion with the usual mean-field elec- 
trostatics and minimized it to derive a modified Poisson- 
Boltzmann equation for potential in the double layer |64| . 
Their model predicts a steady, equilibrium profile of ions 
with a condensed layer at the steric limit for large zeta 
potentials. By extending this approach to obtain the 
chemical potential, Kilic, Bazant and Ajdari have de- 
rived modified-PNP equations and studied steric effects 
on double layer charging, but further model development 
is still needed, not only for steric effects, but also for field- 
dependent permittivity and/or diffusivity. The validity 
of a continuum model at the scale of several atoms in 
the most condensed part of the double layer must also be 
viewed with some skepticism. 

Nevertheless, in spite of these concerns, it is a natural 
first step to study nonlinear electrochemical relaxation 
using the standard PNP equations as we have done. The 
details of our results will surely change with modified 
transport equations and/or surface boundary conditions, 
but we expect that some key features to be robust. For 
example, steric effects could decrease the capacitance of 
the double layer at large zeta potentials, but surface con- 
duction and adsorption, coupled to bulk diffusion, should 
still occur, albeit perhaps with smaller magnitude for a 
given applied field. Also the mathematical aspects of our 
boundary-layer analysis, such as the derivation of the 
surface conservation law l|49(l . could be applied to any 
continuum transport equations. 
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